\section{Code for problem 4}
\label{code4}
\small
\begin{verbatim}
# Vectors for the simulated s(5)'s and q(5)'s
ss <- c()
qs <- c()

# The number of repetitions
N <- 1000
for(j in 1:N)
{
  # The number of steps used in each repetiion
  NrSteps <- 1000
  # The length of the steps
  h <- 5/NrSteps

  # Vectors for the simulated r, s and q
  r <- rep(0,1001)
  s <- rep(0,1001)
  q <- rep(0,1001)
  
  # Initial values for r, s and q.
  r[1] <- 0.044
  s[1] <- 0
  q[1] <- 1
  for(i in 1:NrSteps)
  {
    # We draw a random normal and use this to find r in the next step
    # We also calculate s and q in the next step
    Z <- rnorm(1)
    r[i+1] <- r[i]+0.2*(0.05-r[i])*h+0.1*sqrt(h)*Z+0.01*h*(Z^2-1)/4
    s[i+1] <- s[i]+(1+s[i]*r[i])*h
    q[i+1] <- q[i]+q[i]*r[i]*h
  }

  # We only keep s(5) and q(5)
  ss <- c(ss,s[NrSteps])
  qs <- c(qs,q[NrSteps])
}

# The estimate for the p/q(0) is the mean of the q(5)'s divided by the mean of the s(5)'s
sum(qs)/sum(ss)
# We plot histograms for the s(5)'s and q(5)'s
hist(ss, main="Histogram of simulated E(s(5))'s", xlab="")
hist(qs, main="Histogram of simulated E(q(5))'s", xlab="")
\end{verbatim}

\newpage
